Appendix 3
Katie Florko
11/07/2023
1 Overview
This document contains the code tutorial from “A comparison of statistical models used to characterize species-habitat associations with movement data”. This code tutorial demonstrates the analysis of ringed seal movement in eastern Hudson Bay, Canada.
2 Set-up
2.1 Load packages
# data wrangling
library(here)
library(tidyverse)
library(lubridate)
# mapping
library(rnaturalearth)
library(rgdal)
library(terra)
library(sf)
library(viridis)
# fitting models
library(amt) # for selection functions
library(momentuHMM) # for hidden Markov model2.2 Fish data
Next we load in the fish data. This is a subset of the data from Florko et al. (2021).
# load fish data
fish <- read.csv(here("data/preydiv.csv"))Visualize the fish data.
#prepare world data for mapping
## List of azimuthal equal area - HudsonBay
natearth <- ne_countries(scale = "medium",returnclass = "sp")
natearth <- natearth[which(natearth$continent!="Antarctica"),]
nat_trans <- spTransform(natearth,CRS("+proj=laea +lat_0=60 +lon_0=-85 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))## Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
## to PROJ string## Warning in wkt(obj): CRS object has no comment# plot fish map
fishmap <- ggplot() +
geom_tile(data = fish, aes(x = lon,y = lat,fill = preydiv)) +
scale_fill_viridis(option = "mako", limits = c(0.5, 0.75), name = "Prey\ndiversity") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F)
fishmap
2.3 Seal data
Next we load in the seal movement data. This is a subset of the data from Florko et al. (2023).
# load seal data
seal <- read.csv(here("data/seal_track_m.csv"))
head(seal)## lon lat date id
## 1 357940.9 -368949.6 2012-10-29 14:00 116484_1
## 2 371594.3 -353092.9 2012-10-30 14:00 116484_1
## 3 388429.2 -361478.4 2012-10-31 14:00 116484_1
## 4 399326.2 -357288.1 2012-11-01 14:00 116484_1
## 5 411788.2 -349500.7 2012-11-02 14:00 116484_1
## 6 407708.1 -372564.6 2012-11-03 14:00 116484_1# ensure the data is in the correct format
seal <- seal %>%
mutate(id = as.character(id),
date = as.Date(date)) Visualize the seal data on top of the fish data.
# plot p
sealfishmap <- fishmap +
geom_point(data=seal, aes(x=lon, y=lat), alpha = 0.6, color = "#FCEEAE") +
geom_path(data=seal, aes(x=lon, y=lat), color = "#FCEEAE")
sealfishmap
3 Fit models
We will fit four models: 1) a resource selection function (RSF), 2) a step selection function (SSF), 3) a SSF that allows the habitat covariate to affect movement, and 4) a hidden Markov model (HMM). All four of these models will include prey diversity as a covariate.
3.1 RSF
We will fit the RSF first, as it is the most simple of the four models. In preparation, we will generate an availability sample, create a map to visualize the spatial distribution of the used and available locations, and extract covariates for the used and available locations. Next, we will fit the model and view the model summary.
All of the step selection functions (both the RSF and two SSFs) will be fit using the amt package (Signer et al. 2019).
# prep data and generate availability sample
set.seed(2023)
data_rsf <- seal %>%
make_track(lon, lat, date) %>% # convert data to track format
random_points() # generate availability sample; default is ten times as many available points as observed points
# plot used vs available locations on-top of prey diversity
data_rsf_map <- sealfishmap +
geom_point(data=data_rsf, aes(x=x_, y=y_, color = case_), alpha = 0.6) +
scale_color_manual(values = c("black", "#FCEEAE"),
label = c("Available", "Observed"), name = "Data type")
data_rsf_map
We can see that the availability sample is generated within the minimum convex polygon of the used samples.
# rasterize and extract prey diversity covariate
fish_raster <- terra::rast(fish, crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
data_rsf <- data_rsf %>%
extract_covariates(fish_raster)
# fit rsf (a binomial logistic regression)
rsf1 <- data_rsf %>%
amt::fit_rsf(case_ ~ preydiv, model = TRUE)
# see model summary
summary(rsf1)##
## Call:
## stats::glm(formula = formula, family = stats::binomial(link = "logit"),
## data = data, model = TRUE)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5753 -0.4729 -0.4294 -0.3754 2.3363
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.118 1.400 -4.369 1.25e-05 ***
## preydiv 6.353 2.312 2.748 0.006 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 938.28 on 1539 degrees of freedom
## Residual deviance: 930.60 on 1538 degrees of freedom
## AIC: 934.6
##
## Number of Fisher Scoring iterations: 53.2 SSF
Next we will fit our two SSFs. The workflow is similar to that of the RSF.
# prep data and generate availability sample
set.seed(2023)
data_ssf <- seal %>%
make_track(lon, lat, date) %>% # convert data to track format
steps() %>% # convert track data to step format (i.e., with a start and an end)
random_steps() %>% # generate availability sample
arrange(case_)
# plot used vs available locations on-top of prey diversity
data_ssf_map <- sealfishmap +
geom_point(data=data_ssf, aes(x=x2_, y=y2_, color = case_), alpha = 0.6) +
scale_color_manual(values = c("black", "#FCEEAE"),
label = c("Available", "Observed"), name = "Data type")
data_ssf_map
We can see that the availability sample is generated at each step and is not restricted to the minimum convex polygon.
# extract prey diversity covariate
data_ssf <- data_ssf %>%
extract_covariates(fish_raster, where = "end") #sample at end of step
# transform movement covariates
data_ssf <- data_ssf %>%
mutate(cos_ta = cos(ta_),
log_sl = log(sl_))
# fit ssfs
## ssf1: ssf without covariate affecting movement kernel
ssf1 <- data_ssf %>%
fit_clogit(case_ ~ preydiv + log_sl + cos_ta + strata(step_id_), model = TRUE)
## ssf2: ssf with covariate affecting movement kernel
ssf2 <- data_ssf %>%
fit_clogit(case_ ~ preydiv*log_sl + cos_ta + strata(step_id_), model = TRUE)
# see model summaries
summary(ssf1)## Call:
## coxph(formula = Surv(rep(1, 1518L), case_) ~ preydiv + log_sl +
## cos_ta + strata(step_id_), data = data, model = TRUE, method = "exact")
##
## n= 1516, number of events= 138
## (2 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## preydiv 0.957799 2.605953 6.772711 0.141 0.888
## log_sl -0.008404 0.991631 0.083328 -0.101 0.920
## cos_ta 0.031707 1.032215 0.130643 0.243 0.808
##
## exp(coef) exp(-coef) lower .95 upper .95
## preydiv 2.6060 0.3837 4.477e-06 1.517e+06
## log_sl 0.9916 1.0084 8.422e-01 1.168e+00
## cos_ta 1.0322 0.9688 7.990e-01 1.333e+00
##
## Concordance= 0.494 (se = 0.029 )
## Likelihood ratio test= 0.09 on 3 df, p=1
## Wald test = 0.09 on 3 df, p=1
## Score (logrank) test = 0.09 on 3 df, p=1summary(ssf2)## Call:
## coxph(formula = Surv(rep(1, 1518L), case_) ~ preydiv * log_sl +
## cos_ta + strata(step_id_), data = data, model = TRUE, method = "exact")
##
## n= 1516, number of events= 138
## (2 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## preydiv 2.791e+01 1.318e+12 2.177e+01 1.282 0.200
## log_sl 1.663e+00 5.273e+00 1.286e+00 1.293 0.196
## cos_ta 2.338e-02 1.024e+00 1.311e-01 0.178 0.858
## preydiv:log_sl -2.744e+00 6.433e-02 2.101e+00 -1.306 0.192
##
## exp(coef) exp(-coef) lower .95 upper .95
## preydiv 1.318e+12 7.589e-13 3.919e-07 4.430e+30
## log_sl 5.273e+00 1.897e-01 4.241e-01 6.556e+01
## cos_ta 1.024e+00 9.769e-01 7.917e-01 1.324e+00
## preydiv:log_sl 6.433e-02 1.555e+01 1.046e-03 3.955e+00
##
## Concordance= 0.555 (se = 0.029 )
## Likelihood ratio test= 1.81 on 4 df, p=0.8
## Wald test = 1.81 on 4 df, p=0.8
## Score (logrank) test = 1.81 on 4 df, p=0.83.3 HMM
Finally, we will fit the HMM using momentuHMM (McClintock and Michelot, 2018). In preparation, we define initial parameters, and then update the parameters using our fit model, to ultimately fit a more refined model.
# prepare dataset
data_hmm <- seal %>%
make_track(lon, lat, date) %>%
extract_covariates(fish_raster) %>%
mutate(ID = 1,
x = x_,
y = y_,
date = t_) %>%
dplyr::select(ID, x, y, date, preydiv) %>%
as.data.frame()
# convert into momentuHMM format
data_hmm <- momentuHMM::prepData(data_hmm, coordNames = c("x", "y"),
type = "UTM",
covNames = c("preydiv"))
# define initial parameters
nbStates <- 3 # number of states
stepDist <- "gamma" # step distribution
angleDist <- "vm" # turning angle distribution
mu0 <- c(5000, 18000, 38000) # mean step length for each state
sigma0 <- c(5000, 10000, 10000) # sd step length for each state
stepPar0 <- c(mu0, sigma0)
kappa0 <- c(0.35, 0.55, 0.5) # turning angle for each state
formula = ~ preydiv # identify covariates
# fit basic 3-state hmm
set.seed(2023)
hmm1 <- momentuHMM::fitHMM(data=data_hmm, nbStates=nbStates,
dist=list(step=stepDist,angle=angleDist),
Par0=list(step=stepPar0,angle=kappa0))
# retrieve parameters to refine model
Par0_hmm1 <- momentuHMM::getPar0(hmm1, formula=formula)
# fit a refined hmm with parameters from hhm1
set.seed(2023)
hmm2 <- momentuHMM::fitHMM(data=data_hmm, nbStates=3,
dist=list(step=stepDist,angle=angleDist),
Par0=Par0_hmm1$Par,
delta0 = Par0_hmm1$delta,
beta0 = Par0_hmm1$beta,
formula=formula)We plot the decoded states (estimated behaviours).
# plot decoded states
data_hmm$state <- as.factor(momentuHMM::viterbi(hmm2))
hmmstate_plot <- ggplot() +
scale_fill_viridis(option = "mako", limits = c(0.5, 0.75), name = "Prey\ndiversity") +
geom_path(data=data_hmm, aes(x=x, y=y, color = state, group =ID)) +
geom_point(data=data_hmm, aes(x=x, y=y, color = state, shape = state), size=2, alpha = 0.8) +
scale_color_manual(values = c("#99DDB6", "#539D9C", "#312C66"),
labels = c("Localized\nsearch", "Area-restricted\nsearch", "Travelling"),
name = "HMM state") +
scale_shape_manual(values = c(15, 16, 17),
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling"),
name = "HMM state") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F)
hmmstate_plot
Plot a histogram of step lengths for each state. In order to do this, we need to transform the movement data to UTM and refit our HMM in order to get step length in a meaningful unit (i.e., kilometers)
## prep data in lat/lon and refit model (so step length is in kilometers)
data_hmm <- seal %>%
make_track(lon, lat, date) %>%
extract_covariates(fish_raster) %>%
mutate(ID = 1,
x = x_,
y = y_,
date = t_) %>%
dplyr::select(ID, x, y, date, preydiv) %>%
as.data.frame() %>%
st_as_sf(coords = c("x", "y")) %>%
st_set_crs("+proj=laea +lat_0=60 +lon_0=-85 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0") %>%
st_transform("+proj=longlat +datum=WGS84") %>%
mutate(long = unlist(map(geometry,1)),
lat = unlist(map(geometry,2))) %>%
st_drop_geometry()
# do momentuhmm which calculates step length in KM
data_hmm <- momentuHMM::prepData(data_hmm, coordNames = c("long", "lat"), type = "LL", covNames = c("preydiv"))
# define parameters
nbStates <- 3
stepDist <- "gamma" # step distribution
angleDist <- "vm" # turning angle distribution
mu0 <- c(5, 12, 38)
sigma0 <- c(3, 5, 8)
stepPar0 <- c(mu0, sigma0)
kappa0 <- c(0.35, 0.55, 0.5)
# fit HMM (with step in km)
set.seed(2023)
hmm1_km <- momentuHMM::fitHMM(data=data_hmm, nbStates=nbStates,
dist=list(step=stepDist,angle=angleDist),
Par0=list(step=stepPar0,angle=kappa0))
# get parameters
Par0_hmm1_km <- momentuHMM::getPar0(hmm1_km, formula=formula)
# fit HMM with parameters from hmm1_km
set.seed(2023)
hmm2_km <- momentuHMM::fitHMM(data=data_hmm, nbStates=3,
dist=list(step=stepDist,angle=angleDist),
Par0=Par0_hmm1_km$Par,
delta0 = Par0_hmm1_km$delta,
beta0 = Par0_hmm1_km$beta,
formula=formula)
# add the state estimate from the HMM
data_hmm$state <- as.factor(momentuHMM::viterbi(hmm2_km))
# calculate frequencies of states
v <- momentuHMM::viterbi(hmm2_km)
stateFreq <- table(v) / length(v)
# plot colours
colours.states <- c("#99DDB6", "#539D9C", "#312C66")
# generate sequence for x axis of density functions
x <- seq(0, 50, length=1000)
# get converged mean and sd for each state
meanARS <- hmm2_km$mle$step[1,1]
sdARS <- hmm2_km$mle$step[2,1]
meanCR <- hmm2_km$mle$step[1,2]
sdCR <- hmm2_km$mle$step[2,2]
meanTR <- hmm2_km$mle$step[1,3]
sdTR <- hmm2_km$mle$step[2,3]
# calculate shape and scale of the gamma distributions from mean and sd
sh <- function(mean, sd) { return(mean^2 / sd^2)}
sc <- function(mean, sd) { return(sd^2 / mean)}
# get density functions of the distributions
y_ARS <- dgamma(x, shape=sh(meanARS,sdARS), scale=sc(meanARS,sdARS)) * stateFreq[[1]]
y_CR <- dgamma(x, shape=sh(meanCR,sdCR), scale=sc(meanCR,sdCR)) * stateFreq[[2]]
y_TR <- dgamma(x, shape=sh(meanTR,sdTR), scale=sc(meanTR,sdTR)) * stateFreq[[3]]
# combine densities in a single dataframe for more convenient plotting
df.y_ARS <- data.frame(dens=y_ARS, State="Foraging", x=x)
df.y_CR <- data.frame(dens=y_CR, State="ARS", x=x)
df.y_TR <- data.frame(dens=y_TR, State="Travelling", x=x)
statedis <- rbind(df.y_ARS, df.y_CR, df.y_TR)
# plot distributions
hmm_stepdist_plot <- ggplot() +
geom_line(data=statedis,aes(x=x,y=dens,colour=State,linetype=State), size=1.2) +
scale_colour_manual(values=c(colours.states,"#000000"),
breaks = c('Foraging', 'ARS', 'Travelling'),
labels=c("Localized search", "Area-restricted search", "Travelling")) +
scale_linetype_manual(values=c("solid","solid", "solid"),
breaks = c('Foraging', 'ARS', 'Travelling'),
labels=c("Localized search", "Area-restricted search", "Travelling")) +
ylab("Density") +
xlab("Step length (kms)") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## Warning: Please use `linewidth` instead.hmm_stepdist_plot
Plot a histogram of turning angle for each state.
# generate sequence for x axis of density functions
x <- seq(-pi, pi,length=1000)
# get converged mean and concentration for each state
meanARS <- hmm2_km$mle$angle[1,1]
sdARS <- hmm2_km$mle$angle[2,1]
meanCR <- hmm2_km$mle$angle[1,2]
sdCR <- hmm2_km$mle$angle[2,2]
meanTR <- hmm2_km$mle$angle[1,3]
sdTR <- hmm2_km$mle$angle[2,3]
# get density functions of the distributions
y_ARS <- CircStats::dvm(x, mu=meanARS, kappa=sdARS) * stateFreq[[1]]
y_CR <- CircStats::dvm(x, mu=meanCR, kappa=sdCR) * stateFreq[[2]]
y_TR <- CircStats::dvm(x, mu=meanTR, kappa=sdTR) * stateFreq[[3]]
# combine densities in a single dataframe for more convenient plotting
df.y_ARS <- data.frame(dens=y_ARS, State="Foraging",x=x)
df.y_CR <- data.frame(dens=y_CR, State="ARS",x=x)
df.y_TR <- data.frame(dens=y_TR, State="Travelling",x=x)
cmb <- rbind(df.y_TR, df.y_CR, df.y_ARS)
# plot distributions
hmm_angledist_plot <- ggplot() +
geom_line(data=cmb,aes(x=x,y=dens,colour=State), size = 1.2) +
scale_colour_manual(values=c(colours.states), breaks = c('Foraging', 'ARS', 'Travelling'), labels=c("Localized search", "Area-restricted search", "Travelling")) +
scale_x_continuous(limits=c(-pi,pi))+
ylab("Density") +
xlab("Turn angle (radians) ") +
theme_minimal()
hmm_angledist_plot
4 Prediction line plots
Plotting a model’s estimated relationship between the resource covariate and probability of selection can be useful for general ecological inference. We will calculate the log of the relative selection strength (log-RSS) for each selection function model. The log-RSS is a measure how how likely a step is to end in a proposed location (x1) to one a single reference location (x2), where zero indicates no preference, >1 indicates selection, and <1 indicates avoidance (Avgar et al. 2017, Fieberg et al. 2021).
First, we prepare a dataframe to predict on.
# prep the fish data
newfish <- fish_raster %>%
terra::as.data.frame(xy = TRUE) %>%
filter(x > 100000 & x < 600000 & y > -550000 & y < 0)4.1 RSF
Since the RSF does not incorporate movement, we will calculate the log-RSS of the movement-free habitat selection kernel. This is easily done using log_rss() from amt.
# make a base dataframe to create x1 and x2 from
base <- newfish %>%
mutate(log_sl = log(45),
cos_ta = cos(1))
# x1 is our original dataframe
x1 <- base
# x2 holds prey diversity constant at its mean
x2 <- base %>%
mutate(preydiv = mean(base$preydiv))
# apply log_rss() to each row
log_rss_list <- lapply(1:nrow(x1), function(i) {
# Calculate log-RSS for that row
xx <- log_rss(rsf1, x1[i,], x2[i,], ci = "se")
# Return the element $df
return(xx$df)
})
# combine rows
res1 <- dplyr::bind_rows(log_rss_list)
# plot
line_rsf <- ggplot(res1, aes(x = preydiv_x1, y = (log_rss))) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
geom_line(size = 1, color = "#F8D59F",linetype = 2) +
geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1), alpha = 0.4, fill = "#F8D59F") +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal()
line_rsf
4.2 SSF
We can also calculate the log-RSS of a SSF following the same workflow as for the RSF. Since log_rss() passes to predict(), it is important that the fit SSF includes model = TRUE.
## log-RSS prediction for ssf1
# apply log_rss() to each row
log_rss_list <- lapply(1:nrow(x1), function(i) {
# Calculate log-RSS for that row
xx <- log_rss(ssf1, x1[i,], x2[i,], ci = "se")
# Return the element $df
return(xx$df)
})
# combine rows
res2 <- dplyr::bind_rows(log_rss_list) %>%
mutate(Speed = "without int.")
# plot
line_ssf1 <- ggplot() +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
geom_line(data = res2, aes(x = preydiv_x1, y = log_rss, color = Speed, group = Speed), size = 1, linetype = 3, color = "black") +
geom_ribbon(data = res2, aes(ymin=lwr, ymax=upr, x=preydiv_x1, fill = Speed, group = Speed), alpha = 0.3, fill = "black") +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal()
line_ssf1
For ssf2 (the model where prey diversity interacts with movement), we will estimate the log-RSS for three different step-lengths (slow, moderate, fast). We set these speeds as the 25th, 50th, and 75th percentiles of step-length, then loop the log-RSS for each speed.
## log-RSS prediction for ssf2
# determine the 25th (slow), 50th (moderate), and 75th (fast) percentiles of step-length
nums = seal %>%
make_track(lon, lat, date) %>%
steps %>%
mutate(log_sl = log(sl_)) %>%
summarize(quants = quantile(log_sl, c(0.25, 0.5, 0.75))) %>%
pull()
# set-up to run function for each speed
results_ssf2 <- lapply(nums, function(i) {
x1$log_sl <- i
# calculate log-RSS
log_rss_list <- lapply(1:nrow(x1), function(i) {
# Calculate log-RSS for that row
xx <- log_rss(ssf2, x1[i,], x2[i,], ci = "se")
# Return the element $df
return(xx$df)
})
# bind rows within each speed's prediction
res3 <- dplyr::bind_rows(log_rss_list)
} )
# bind rows of all speed's prediction
results_ssf2 <- dplyr::bind_rows(results_ssf2) %>%
mutate(log_sl_x1 = as.factor(round(log_sl_x1,1)),
Speed = dplyr::case_when(as.factor(log_sl_x1) == '8.4' ~ "slow",
as.factor(log_sl_x1) == "9.2" ~ "moderate",
as.factor(log_sl_x1) == "9.6" ~ "fast"))
# plot
line_ssf2 <- ggplot(results_ssf2, aes(x = preydiv_x1, y = (log_rss))) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
geom_line(size = 1, aes(color = Speed, group = Speed, linetype = Speed)) +
geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1, fill = Speed, group = Speed), alpha = 0.3) +
scale_colour_manual(values=c(colours.states,"#000000")) +
scale_fill_manual(values=c(colours.states,"#000000")) +
scale_linetype_manual(values = c("solid", "solid", "solid")) +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal()
line_ssf2
Typically these models would be interpreted independently. But it’s worth noting that while the effects are minimal and the confidence intervals overlap, when comparing ssf1 and ssf2, the ssf2 provides more information about the animal’s relationship with prey diversity.
# Plot
ggplot(results_ssf2, aes(x = preydiv_x1, y = (log_rss))) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
geom_line(size = 1, aes(color = Speed, group = Speed, linetype = Speed)) +
geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1, fill = Speed, group = Speed), alpha = 0.3) +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal() +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
geom_line(data = res2, aes(x = preydiv_x1, y = log_rss, color = Speed, group = Speed, linetype = Speed), size = 1) +
geom_ribbon(data = res2, aes(ymin=lwr, ymax=upr, x=preydiv_x1, fill = Speed, group = Speed), alpha = 0.2) +
scale_color_manual(values = c("#99DDB6", "#539D9C", "#312C66", "black"),
labels = c("Fast", "Moderate", "Slow", "Not consid."),
name = "Speed") +
scale_fill_manual(values = c("#99DDB6", "#539D9C", "#312C66", "black"),
labels = c("Fast", "Moderate", "Slow", "Not consid."),
name = "Speed") +
scale_linetype_manual(values = c("solid", "solid", "solid", "dotted"),
labels = c("Fast", "Moderate", "Slow", "Not consid."),
name = "Speed") 
4.3 HMM
Finally, we grab and plot the stationary state probabilities. This is easily done using plotStationary() from momentuHMM.
# grab stationary probabilities
ps <- momentuHMM::plotStationary(hmm2, plotCI= TRUE, return = TRUE)state1 <- ps$preydiv$'state 1' %>% mutate(state = 1)
state2 <- ps$preydiv$'state 2' %>% mutate(state = 2)
state3 <- ps$preydiv$'state 3' %>% mutate(state = 3)
# bind to one data frame
pdat <- rbind(state1, state2, state3) %>%
mutate(state = as.character(state))# plot
line_hmm <- ggplot() +
geom_line(data = pdat, aes(x = cov, y = est, color = state)) +
geom_ribbon(data = pdat, aes(x=cov, y=est, ymax=est+se, ymin=est-se, fill = state),
alpha = 0.4, show.legend = TRUE) +
ylab("Stationary state probabilties") +
xlab("Prey diversity") +
scale_color_manual(values = c("#99DDB6", "#539D9C", "#312C66"), name = "HMM state",
labels=c("Localized search", "Area-restricted search", "Travel")) +
scale_fill_manual(values = c("#99DDB6", "#539D9C", "#312C66"), name = "HMM state",
labels=c("Localized search", "Area-restricted search", "Travel")) +
theme_minimal()
line_hmm5 Prediction maps
Now we will estimate the utilization distribution from each model to demonstrate the inference from each model. The utilization distribution is defined as the two-dimensional relative frequency distribution of space use of an animal (Van Winkle 1975). This is a simple calculation for the RSF, where we multiply the model coefficient with the resource (prey diversity), exponentiate (since it is a logistic regression), and normalize the estimate. The calculations are more complex for the SSFs since they are conditional models that integrate the movement process. Thus, for the SSFs we calculate the steady-state utilization distribution (SSUD), which is the long-term expectation of the space-use distribution across the landscape (Signer et al. 2017). amt has functions to estimate the SSUD.
5.1 RSF
We can predict the estimated probability of use from the RSF by hand.
# grab model coefficients
modcoef <- summary(rsf1)$coef
# prediction for each cell
x <- exp(modcoef[1] + (modcoef[2] * newfish$preydiv))
# scale the data
x <- x / (1 + x)
# range fn
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
# set the range from zero to one
newfish$rsf_prediction <-range01(x)
# plot
map_rsf <- ggplot() +
geom_tile(data = newfish, aes(x = x,y = y, fill = rsf_prediction)) +
scale_fill_viridis(option = "mako", name = "RSF prediction", limits = c(0, 0.3)) +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) +
theme_void()## Regions defined for each Polygonsmap_rsf
5.2 SSF
We can use the amt functions to estimate the SSUDs from the simple SSF that does not allow prey diversity to affect the movement kernel.
# generate availability sample
set.seed(2023)
data_ssf <- seal %>%
mutate(date = as.POSIXct(date)) %>%
make_track(lon, lat, date) %>%
steps() %>%
random_steps() %>%
arrange(case_) %>%
amt::extract_covariates(fish_raster, where = "both") # %>%
#na.omit()
# fit SSF1 model
m1 <- data_ssf |>
fit_clogit(case_ ~ preydiv_end + cos(ta_) + log(sl_) +
strata(step_id_))We will now simulate a track and visually observe it.
# set starting position for the simulation
start <- make_start((seal %>%
mutate(date = as.POSIXct(date)) %>%
make_track(lon, lat, date))[1,])
# Set constants
n_steps = 1e3 # number of steps
n_steps1 = n_steps + 1 # number of steps +1
burnin <- n_steps/50 # number of locations to remove for burn-in
# generate redistribution kernel
k1 <- redistribution_kernel(m1, map = fish_raster, start = start,
stochastic = TRUE,
tolerance.outside = 1,
n.control = 1e3)
# simulate path
set.seed(2023)
p1 <- amt::simulate_path(x = k1, n = n_steps, start = start, verbose = TRUE)
# burn-in
p1_burnt <- p1 %>% slice(-c(1:burnin))
# plot simulated track
ssf_track_1 <- fishmap +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
geom_point(data = p1_burnt, aes(x = x_,y = y_), alpha = 0.61) +
geom_path(data = p1_burnt, aes(x = x_,y = y_)) +
theme_void()## Regions defined for each Polygonsssf_track_1
We can see that it mostly stays within the study area. We will use this track to estimate the SSUD, and visualize the results.
# estimate SSUD
uds_ssf1 <- tibble(rep = 1:n_steps1,
x_ = p1$x_, y_ = p1$y_,
t_ = p1$t_, dt = p1$dt) |>
filter(!is.na(x_)) |>
make_track(x_, y_) |>
hr_kde(trast = fish_raster, which_min = "global") %>%
hr_ud() %>%
terra::as.data.frame(xy = TRUE)
# plot SSUD
map_ssf1 <- ggplot() +
geom_tile(data = uds_ssf1, aes(x = x,y = y, fill = preydiv)) +
scale_fill_viridis(option = "mako", name = "SSF1 prediction") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) +
theme_void()## Regions defined for each Polygonsmap_ssf1
We will follow the same steps to generate a SSUD for the SSF that allows prey diversity to affect the movement kernel.
# fit SSF2 model
m2 <- data_ssf |>
fit_clogit(case_ ~ preydiv_end + cos(ta_)+
preydiv_end:log(sl_) +
strata(step_id_))
# set starting position for the simulation
set.seed(2023)
start <- make_start((seal %>%
mutate(date = as.POSIXct(date)) %>%
make_track(lon, lat, date))[1,])
# generate redistribution kernel
k2 <- redistribution_kernel(m2, map = fish_raster, start = start,
stochastic = TRUE,
tolerance.outside = 1,
n.control = 1e3)
# Now simulate a path of length 1e3
n_steps = 1e3
n_steps1 = n_steps + 1
set.seed(2023)
p2 <- amt::simulate_path(x = k2, n = n_steps, start = start)
# plot simulated track
ssf_track_2 <- fishmap +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
geom_point(data = p2, aes(x = x_,y = y_), alpha = 0.61) +
geom_path(data = p2, aes(x = x_,y = y_)) +
theme_void()## Regions defined for each Polygonsssf_track_2
# estimate SSUD
uds_ssf2 <- tibble(rep = 1:n_steps1,
x_ = p2$x_, y_ = p2$y_,
t_ = p2$t_, dt = p2$dt) |>
filter(!is.na(x_)) |>
make_track(x_, y_) |>
hr_kde(trast = fish_raster, which_min = "local") %>%
hr_ud() %>%
terra::as.data.frame(xy = TRUE)
# plot SSUD
map_ss2 <- ggplot() +
geom_tile(data = uds_ssf2, aes(x = x,y = y, fill = preydiv)) +
scale_fill_viridis(option = "mako", name = "SSF2 prediction") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) +
theme_void()## Regions defined for each Polygonsmap_ss2
5.3 HMM
We will first estimate the stationary state probabilities of each state based on prey diversity. This is easily done using the momentuHMM function stationary().
# grab estimated stationary state probabilities from our fitted HMM
x <- as.data.frame(momentuHMM::stationary(hmm2, data.frame(preydiv = newfish$preydiv)))
newfish$hmm_state1 <- x$state.1
newfish$hmm_state2 <- x$state.2
newfish$hmm_state3 <- x$state.3
# prepare data
newfish_long <- newfish %>%
tidyr::pivot_longer(cols = hmm_state1:hmm_state3,
names_to = "model", values_to = "prediction") %>%
mutate(dplyr::across(model, factor, levels=
c("hmm_state1", "hmm_state2", "hmm_state3")))
# plot
map_hmm <- ggplot() +
geom_tile(data = newfish_long, aes(x = x, y = y, fill = prediction)) +
scale_fill_viridis(option = "mako", limits = c(0,1)) +
labs(fill = 'HMM predicted\nprobability') +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), colour = "white", fill = "grey90", size = 0.3) +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) +
facet_wrap(~model, labeller = as_labeller(c('hmm_state1' = "Localized search",
'hmm_state2' = "Area-restricted search",
'hmm_state3' = "Travel"))) +
theme_void()## Regions defined for each Polygonsmap_hmm
6 References
Avgar, T., Lele, S. R., Keim, J. L., & Boyce, M. S. (2017). Relative selection strength: Quantifying effect size in habitat‐and step‐selection inference. Ecology and Evolution, 7(14). 5322-5330.https://doi.org/10.1002/ece3.3122
Fieberg, J., Signer, J., Smith, B., & Avgar, T. (2021). A ‘How to’guide for interpreting parameters in habitat‐selection analyses. Journal of Animal Ecology, 90(5), 1027-1043. https://doi.org/10.1111/1365-2656.13441
Florko, K.R.N., Tai, T.C., Cheung, W.W.L., Sumaila, U.R., Ferguson, S.H., Yurkowski, D.J., Auger-Méthé, M. (2021). Predicting how climate change threatens the prey base of Arctic marine predators. Ecology Letters, 24: 2563-2575. https://doi.org/10.1111/ele.13866
Florko, K.R.N., Shuert, C.R., Cheung, W.W.L., Ferguson, S.H., Jonsen, I.D., Rosen, D.A.S., Sumaila, U.R., Tai, T.C., Yurkowski, D.J., Auger-Méthé, M. (2023). Linking movement and dive data to prey distribution models: new insights in foraging behavior and potential pitfalls of movement analyses. Movement Ecology, 11:17 https://doi.org/10.1186/s40462-023-00377-2
McClintock, B.T., Michelot, T. (2018). momentuHMM: R package for generalized hidden markov models of animal movement. Methods in Ecology and Evolution, 9, 1518-1530. https://doi.org/10.1111/2041-210X.12995
Signer, J., Fieberg, J., & Avgar, T. (2017). Estimating utilization distributions from fitted step‐selection functions. Ecosphere, 8(4), e01771. https://doi.org/10.1002/ecs2.1771
Signer, J., J. Fieberg, and T. Avgar. (2019). Animal movement tools (amt): R package for managing tracking data and conducting habitat selection analyses. Ecology and Evolution, 9:880–890. https://doi.org/10.1002/ece3.4823
Van Winkle, W. (1975). Comparison of several probabilistic home-range models. The Journal of wildlife Management, 118-123. https://doi.org/10.2307/3800474